
### CODE FOR VARIATIONAL NETWORK MODEL TEST
##
## From Barbera's original replication archive for the US test, we have the following files:
## "adj-matrix-US-rdata","results-elites-US-rdata","elites-data.Rdata", and "users-data-US.rdata"
##
##


### ########
### Packages
### ########

setwd("~/Documents/working/network/simulation")
rm(list = ls())
library("emIRT")
library("Matrix")
library("methods")

### #####################
### Common Est Parameters
### #####################

#cMaxIter <- 5e3
cMaxIter <- 1500
cThresh <- 1e-3
cUsers <- 10000

### ############
### Load US Data
### ############

load("adj-matrix-US.rdata")
load("results-elites-US.rdata")
country <- 'US'
load("elites-data.Rdata")

# US:
if (country=="US"){
	us <- elites.data[['US']]
	parties <- merge(
		data.frame(screen_name = colnames(y), stringsAsFactors=F),
		us[,c("screen_name", "party")], sort=FALSE, all.x=TRUE)$party
	start.phi <- rep(0, length(parties))
	start.phi[parties=="D"] <- -1
	start.phi[parties=="R"] <- 1
}


### ##################################
### US -- Convert from Sparse to Dense
### ##################################

inform <- which(rowSums(y)>10)
set.seed(12345)
subset.i <- sample(inform, cUsers)
y2 <- y[subset.i, ]
start.phi <- start.phi[which(colSums(y2)>200)]
y2 <- y2[,which(colSums(y2)>200)]
y2 <- as.matrix(y2)

### ##################
### Identify US Anchor
### ##################

results.subset <- results[which(colSums(y2)>200),]
idxAnchor <- which(results.subset$phi == max(results.subset$phi))[1]

### ##################
### US -- Get Key Qtys
### ##################

cN <- nrow(y2)
cJ <- ncol(y2)

lP <- makePriorsE(cN,
                  cJ
                  )

set.seed(1)

lS <- getStartsE(cN,
                 cJ,
                 y2
                 )


time <- system.time({
lEus <- networkIRT(.y = y2,
                   .starts = lS,
                   .priors = lP,
                   .control = {list(verbose = TRUE,
                                    maxit = cMaxIter,
                                    convtype = 2,
                                    thresh = cThresh,
                                    threads = 8
                                    )
                           },
                   .anchor_item = idxAnchor
                   )
})

time[3]  #2094 seconds

### ##########
### US Results
### ##########
dfVI <- data.frame(screen_name = colnames(y2),
                   location = lEus$means$w,
                   popularity = lEus$means$alpha
                   )

dfVI2 <- merge(dfVI,
               results.subset,
               by = "screen_name",
               suffixes = c("", ".pb")
               )

cor(dfVI2$phi, dfVI2$location)

## Save intermediate output
#save(lEus, time, dfVI, dfVI2, file="US_vi_out.rda")
#load("US_vi_out.rda")

### Graphics
load("users-data-US.rdata")
matches <- match(users$uid, rownames(y2))
users$theta.vi <- lEus$means$theta[matches]
users <- users[!is.na(matches),]

dfVI2$phi <- (dfVI2$phi - mean(dfVI2$phi,na.rm=TRUE))/sd(dfVI2$phi,na.rm=TRUE)
dfVI2$location <- -1*(dfVI2$location - mean(dfVI2$location,na.rm=TRUE))/sd(dfVI2$location,na.rm=TRUE)
users$theta <- (users$theta - mean(users$theta))/sd(users$theta)
users$theta.vi <- -1*(users$theta.vi - mean(users$theta.vi))/sd(users$theta.vi)

cor(users$theta, users$theta.vi, use="pairwise")
cor(dfVI2$phi, dfVI2$location, use="pairwise")

pdf("elites_barbera.pdf")
par(mar=c(5, 5, 4, 2) + 0.1)
plot(dfVI2$phi,
     dfVI2$location,
     xlab = "Barbera (2015)",
     ylab = "EM Estimate",
     main = "Elite Positions", cex.main=1.5, cex.lab=1.7, cex.axis=1.5, ylim=c(-2,2), type="n")
text(1,-1,expression(paste(rho," = 0.99")), cex=2)
abline(a=0, b=1, lwd=4, col="grey")
points(dfVI2$phi, dfVI2$location, pch=16, cex=1.2)
dev.off()

pdf("voter_barbera.pdf")
par(mar=c(5, 5, 4, 2) + 0.1)
plot(users$theta,
     users$theta.vi,
     xlab = "Barbera (2015)",
     ylab = "EM Estimate",
     main = "Voter Positions", cex.main=1.5, cex.lab=1.7, cex.axis=1.5, type="n")
text(1.5,-2, expression(paste(rho," = 0.98")), cex=2)
abline(a=0, b=1, lwd=4, col="grey")
points(users$theta, users$theta.vi, pch=16, cex=0.7)
dev.off()




